function sT0(f,fi)
 global Ti T0i T0 wi w wi0 C2
   
   fxy=ave(convect(Ti));
%    fxy = C2*fxy + 1e-6*T0i+0.05*dif3(T0i);
   fxy = fxy - dif3(T0i);
%    fxy = C2*fxy - Lpara(T0i);

   work=T0*f+T0i*fi-fxy;
  %%
   if(f<0.6)%the 1st time
       T0=T0i;
       T0i=work;
   else %the 2nd time
       T0i=work;
   end
end